計算機で一様擬似乱数を生成する．
実数 \(a < b\) に対し，区間 \([a, b]\) 上の 一様分布 (uniform distribution) とは，確率密度関数 \[p(x) = \begin{cases} 0 & \text{if} \ \ x < a \\ \frac{1}{b - a} & \text{if} \ \ a \leq x \leq b \\ 0 & \text{if} \ \ x > b \end{cases}\tag{2.1}\] を持つ確率分布．\(\mathcal{U}[a, b]\) と書く．区間を指定せずに一様分布というと，\(\mathcal{U}[0, 1]\) を指す．
set.seed(1) # seedの設定
runif(3) # 長さ3の一様擬似乱数## [1] 0.2655087 0.3721239 0.5728534関数 runif の出力は一様乱数と認識しても問題ないか？
\(\rightsquigarrow\) スペクトル (Spectral) 検定，コルモゴロフ・スミルノフ (Kolmogorov-Smirnov) 検定 を適用．
# スペクトル検定
set.seed(1)
x <- runif(1000)
y <- runif(1000)
z <- runif(1000)
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)立方体に均等に分布しているように見える．
# コルモゴロフ・スミルノフ検定
set.seed(1)
u <- runif(1000)
ks.test(u, punif)##
## One-sample Kolmogorov-Smirnov test
##
## data: u
## D = 0.024366, p-value = 0.5928
## alternative hypothesis: two-sided一様分布から乖離しているのであれば，\(p\)-値は小さくなるはずであるが，ここでは一様擬似乱数を観測として見たときの \(p\)-値は小さくない．
\(\rightsquigarrow\) 一様擬似乱数 (初期設定はMersenne-Twister法) は一様分布からの独立な乱数列と見ても顕著な差が見られない．
\(a, b, y_0, n\): 事前に決める整数．
\[ y_m = (a y_{m - 1} + b) \mod n \ (m = 1, 2, \ldots) \\ x_m = \frac{y_m}{n}\] として\(\{ x_m: m = 1, 2, \ldots \}\) を出力する．
u <- numeric(3*1e3)
y <- 1234
n <- 2^31 - 1
a <- 65539
b <- 0
for(i in 1:length(u)){
y <- (a*y + b) %% n
u[i] <- y/n
}
u_mat <- matrix(u, nrow=3)
x <- u_mat[1,]
y <- u_mat[2,]
z <- u_mat[3,]
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)2次元平面に整列している様子が見て取れる．
以下，理論的な記述の際は乱数を扱い，実際の計算では擬似乱数を用いるが，擬似乱数も乱数と表記する．
ここで，\(\mathcal{U}[0, 1]\) からの一様乱数が生成できるとする．\(X \sim \mathcal{U}[0, 1]\) に対して，線形変換
\[ Y = b + (a - b) X \] を施すことで，\(\mathcal{U}[a, b]\) に従う乱数が生成できる．
\(\because\) \(\mathcal{U}[a, b]\) は累積分布関数 \[ F(x) = \frac{x - a}{b - a} \ (a \leq x \leq b) \] を持つが，線形変換してできた \(Y\) の従う分布も \[ \mathbb{P}(Y \leq x) = \mathbb{P}\left( X \leq \frac{x - b}{a - b} \right) = \frac{x - b}{a - b} = F(x) \ (a \leq x \leq b) \] となり，同じ累積分布関数を持つため．
一様乱数をもとに，与えられた1次元の確率分布 \(P\) に従う乱数を生成する．
累積分布関数を \(F\) とすると， \[ \mathbb{P}(X \leq x) = F(x) \ (-\infty < x < \infty) \] となるような \(X\) を生成すればよい．
連続な確率分布であれば，任意の \(u \in (0, 1)\) に対して，\(F(x) = u\) となる \(x\) がただ1つに決まる．
\(\rightsquigarrow\) \(F^{-1}(u) = x\) で逆関数 \(F^{-1}\) が定義できる．このとき，一様乱数 \(U\) により， \[ X = F^{-1}(U) \] で \(P\) に従う確率変数を作れる．
\[\because \mathbb{P}(X \leq x) = \mathbb{P}(F^{-1}(U) \leq x) = \mathbb{P}(U \leq F(x)) = F(x) \] このようにして確率分布 \(P\) に従う乱数を生成する方法を 逆変換法 (inversion method) という．
\(\lambda > 0\) とする．指数分布 \(\mathcal{E}(\lambda)\) の累積分布関数は， \[ F(x) = 1 - e^{-\lambda x} \] である．したがって，一様乱数 \(U\) を用いて \[ F^{-1}(U) = - \frac{\log (1 - U)}{\lambda} \] によって指数乱数が発生できる．ここで，\(1 - U \sim \mathcal{U}[0, 1]\) より，\(-\log(U)/\lambda\) によっても指数乱数を生成できる．
R 言語で実装する．パラメータ \(\lambda\) を \(\lambda = 2.0\) とする．
set.seed(1)
lambda <- 2.0
u <- runif(3)
-log(u)/lambda## [1] 0.6630539 0.4942642 0.2785628組み込み関数 rexp でも生成する．
set.seed(1)
lambda <- 2.0
rexp(3, rate=lambda)## [1] 0.37759092 0.59082139 0.07285336同じ seed を用いても異なる値が出力されていることから，異なった実装がなされているようである．
計算時間の比較を行う．
lambda <- 2.0
n <- 1e5
microbenchmark(A <- -log(runif(n))/lambda, times=100)## Unit: milliseconds
## expr min lq mean median uq
## A <- -log(runif(n))/lambda 2.932916 3.165062 3.313722 3.250917 3.386003
## max neval
## 5.024328 100microbenchmark(B <- rexp(n, rate=lambda), times=100)## Unit: milliseconds
## expr min lq mean median uq
## B <- rexp(n, rate = lambda) 3.465591 3.550562 3.755103 3.663531 3.852399
## max neval
## 5.572789 100実行時間に大きな差はないようである．
実数 \(\mu\) と，正の実数 \(\sigma\) に対し，コーシー分布 (Cauchy distribution) \(\mathcal{C}(\mu, \sigma)\) は，確率密度関数 \[ p(x; \mu, \sigma) = \frac{1}{\pi \sigma \left( 1 + \frac{(x - \mu)^2}{\sigma^2} \right)} \ (-\infty < x < \infty) \] を持つ．\(\mu = 0, \sigma = 1\) としたコーシー分布 \(\mathcal{C}(0, 1)\) に注目する．累積分布関数は， \[ F(x) = \int_{-\infty}^x \frac{1}{\pi (1 + y^2)}dy \] となる．学部の解析学を思い出すと，この積分は， \[F(x) = \frac{\mathrm{arctan} x}{\pi} + \frac{1}{2}\] となる．よって， \[F^{-1}(u) = \tan \left( \pi u - \frac{\pi}{2} \right)\] を得る．したがって，一様乱数 \(U\) を用いて， \[X = \tan (\pi U - \pi/2)\] でコーシー乱数を生成できる．
また，コーシー分布は独立な2つの標準正規乱数 \(X_1, X_2\) の比 \[\frac{X_2}{X_1} \tag{2.2}\] の分布と等しい．以下で，逆変換法，(2.2)の方法，組み込み関数の3つの計算時間を比較する．
n <- 1e5
microbenchmark(A <- tan(pi*runif(n)-pi/2), times=100)## Unit: milliseconds
## expr min lq mean median uq
## A <- tan(pi * runif(n) - pi/2) 4.955723 5.177404 5.373381 5.303579 5.505346
## max neval
## 6.806326 100microbenchmark(B <- rnorm(n)/rnorm(n), times=100)## Unit: milliseconds
## expr min lq mean median uq max
## B <- rnorm(n)/rnorm(n) 10.00244 10.32055 10.66322 10.46786 10.88738 12.83597
## neval
## 100microbenchmark(C <- rcauchy(n), times=100)## Unit: milliseconds
## expr min lq mean median uq max neval
## C <- rcauchy(n) 5.188717 5.419534 5.627272 5.552008 5.713074 7.21623 100上述の逆変換法では，累積分布関数 \(F\) の逆関数 \(F^{-1}\) が存在することを仮定した．
離散の場合は逆関数を持たないが，一般に累積分布関数 \(F(x) = \mathbb{P}((-\infty, x])\) は次の性質を持つことを思い出す．
逆に，上記の性質を満たす関数 \(F\) は \(\mathbb{R}\) 上のある確率分布の累積分布関数である．
ここで，関数 \(F\) が単調非減少で右連続のとき， \[F^{-}(u) = \inf \{x \in \mathbb{R}: u \leq F(x) \}\] によって広義の逆関数 \(F^{-}: \mathbb{R} \to \mathbb{R}\) を定める．広義の逆関数 \(F^{-}\) に対しては，\(F^{-}(F(x)) = x\) や \(F(F^{-}(u)) = u\) は必ずしも成り立たないが，次が成り立つ．
\(x \in \mathbb{R}, u \in [0, 1]\) に対して， \[F^{-}(F(x)) \leq x\] \[ F(F^{-}(u)) \geq u \]
確率分布 \(P\) の累積分布関数 \(F\) の広義の逆関数 \(F^{-}\) を用いることで，連続の場合と同様に \(U \sim \mathcal{U}[0, 1]\) から \(F^{-}(U)\) を計算することで，\(P\) に従う乱数を生成できる．
有限集合 \(-\infty < x_1 < \cdots < x_n < \infty\) にのみ値をとる確率分布を考える．この分布の累積分布関数 \(F\) の広義の逆関数 \(F^{-}\) は， \[ F^{-}(u) = \begin{cases} x_1 & \text{if} \ \ 0 \leq u \leq F(x_1) \\ x_2 & \text{if} \ \ F(x_1) < u \leq F(x_2) \\ \vdots \\ x_n & \text{if} \ \ F(x_{n - 1}) < u \leq 1 \end{cases} \] である．可算無限集合でも同様．